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Abstract:  Domain  decomposition  is  a  class  of  techniques  that  are  designed  to  solve  elliptic  problems 
on  irregular  domains  and  on  multiprocessor  systems.  Typically,  a  domain  is  decomposed  into  many 
smaller  regular  subdomains  and  the  capacitance  system  governing  the  interface  unknowns  is  solved 
by  some  version  of  the  preconditioned  conjugate  gradient  method.  In  this  paper,  we  show  that  for  a 
simple  model  problem  —  Poisson’s  equation  on  a  rectangle  decomposed  into  two  smaller  rectangles 
—  the  capacitance  system  can  be  inverted  exactly  by  fast  Fourier  transform.  No  iteration  is  needed. 
An  exact  eigen-decomposition  of  the  capacitance  matrix  also  makes  possible  a  comparison  of  various 
preconditioners  that  have  been  proposed  in  the  literature.  For  example,  we  show  that  in  the  limit 
as  the  aspect  ratio  of  the  two  rectangles  tend  to  infinity,  the  preconditioner  proposed  by  Golub- 
Mayers  becomes  exact,  but  the  one  proposed  by  Dryja  does  not.  Both  preconditioners,  however, 
are  poor  when  the  aspect  ratio  is  small.  f  _ 
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1.  Introduction 

Domain  decomposition  is  a  class  of  techniques  that  is  designed  to  solve  elliptic  problems  on 
irregular  domains  and  on  multiprocessor  systems.  Typically,  a  domain  is  decomposed  into  many 
smaller  regular  subdomains  and  the  capacitance  system  governing  the  interface  unknowns  is  solved. 
This  is  a  relatively  old  idea  and  can  be  traced  to  the  Schwarz  alternating  procedure  [9], 

Such  methods  are  attractive  in  many  situations.  In  fact,  the  main  reason  for  the  resurgence  of 
this  old  idea  is  its  obvious  advantage  in  implementation  on  multiprocessor  systems.  Even  in  a  se¬ 
quential  computing  environment,  a  natural  partition  of  the  computational  domain  often  exists,  such 
as  in  dividing  a  domain  with  irregular  geometry  into  regular  subregions  for  which  fast  solvers  exist, 
or  in  dividing  a  problem  with  discontinous  coefficients  into  subregions  with  constant  coefficients. 
For  this  and  other  reasons,  domain  decomposition  has  received  a  lot  of  interest  recently. 

Since  the  capacitance  system  is  expensive  to  evaluate  and  expensive  to  solve  by  direct  methods, 
many  of  the  methods  proposed  so  far  in  the  literature  employ  some  form  of  preconditioned  conjugate 
gradient  (PCG)  method  for  its  solution.  In  each  iteration,  the  product  of  the  capacitance  matrix 
and  a  given  vector  is  required,  which  can  be  evaluated  by  solving  problems  on  the  subdomains. 
To  minimize  the  number  of  subdomain  solves,  it  is  imperative  to  have  a  good  preconditioner  for 
the  capacitance  matrix.  Dryja[5]  is  among  the  first  to  introduce  such  a  preconditioner  for  two 
dimensional  problems,  which  is  in  the  form  of  a  pseudo  differential  operator,  namely,  the  square 
root  of  the  one  dimensional  discrete  Laplace  operator.  Later,  Golub  and  Mayers[8]  proposed  a 
modification  which  significantly  reduces  the  number  of  PCG  iterations  needed.  Many  other  methods 
have  been  proposed  along  this  approach,  among  which  we  mention  [1,  2,  3,  6,  7,  8]. 

In  this  paper,  we  show  that  for  a  simple  model  problem  —  Poisson’s  equation  on  a  rectangle 
decomposed  into  two  smaller  rectangles  —  the  capacitance  system  can  be  inverted  exactly  by  fast 
Fourier  transform.  No  PCG  iteration  ia  needed.  We  derive  an  exact  eigen-decomposition  of  the 
capacitance  matrix  which  makes  possible  a  comparison  of  various  preconditioners  that  have  been 
proposed  in  the  literature.  For  example,  we  show  that  in  the  limit  as  the  aspect  ratio  of  the  two 
rectangles  tend  to  infinity,  the  preconditioner  proposed  by  Golub-Mayers  becomes  exact,  but  the 
one  proposed  by  Dryja  does  not.  Both  preconditioners,  however,  are  poor  when  the  aspect  ratio  is 
small. 

In  Section  2,  we  introduce  the  model  problem  and  derive  the  capacitance  system  for  the 
interface.  The  eigen-decomposition  of  the  capacitance  matrix  is  derived  in  Section  3.  Comparisons 
of  various  preconditioners  are  discussed  in  Section  4  and  we  close  in  Section  5  with  some  remarks 
about  extensions  to  irregular  regions  and  divisions  into  more  subdomains. 
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2.  Ths  Modal  Problam  and  the  Interface  System 
Consider  the  following  Poisson’s  equation  : 


A  u  =  /  on  fl 


with  boundary  condition 


v  —  g  on  dO 


and  where  the  domain  fl  is  as  illustrated  in  Fig.  1. 


y 


(2.1) 


Fig.  1  The  domain  Q  and  its  partition 

We  partition  fl  into  two  subdomains  fit  and  fl],  with  a  common  interface  03.  We  use  a  uniform 
mesh  with  grid  size  h  on  fl  with  n  internal  grid  points  in  the  z-direction,  i.e., 

h=  .  1  .. 

(n+  1) 

We  assume  that  /1  and  I2  are  integral  multiples  of  h,  with  mi  internal  grid  points  in  fit  in  the 
y-direction  and  m2  internal  grid  points  in  fl2,  i.e., 

11  —  (mi  +  1  )h 

12  —  (mj  +  1  )h  . 

Consider  a  standard  S-point  centered  difference  approximation  to  (2.1).  If  we  order  the  unknowns 
in  fli  first,  then  those  in  fl2  and  finally  those  in  fl3,  then  the  discrete  solution  vector  u  =  (uj,  U2,  U3) 


satisfies  the  following  linear  system 

(M\ 

0 

a,3\  ( 

^  /  /1 A 

0 

A22 

A23 

“21  =  1/2 

(2.2) 

Ufa 

•*23 

•*33  /  V 

“3  /  \h) 

where  the  matrices  An,  A22,  and  A33  correspond  to  the  discrete  Laplacian  on  fli,  fl2  and  fl3.  and 
A13  and  A23  correspond  to  the  coupling  between  the  unknowns  in  fli  and  fl]  with  those  in  fl3- 


Applying  block  Gaussian  Elimination  to  (2.2),  we  obtain  the  following  system  for  the  interface 
unknowns  U3  : 

CU3  =  f 3  ~  A13A11  f\  —  AjjjAjj1  fi  (2-3) 

where 

C  s  A33  —  AijjAjj1  A13  —  A23A22  A23.  (2-4) 

The  right  hand  side  of  (2.3)  can  be  evaluated  by  solving  two  subdomain  problems,  one  each  on  Hi 
and  flj. 

Note  that  C  is  expensive  to  evaluate  because  it  requires  2n  subdomain  solves.  Moreover,  it  is 
generally  dense  and  therefore  a  direct  method  for  solving  (2.3)  could  be  prohibitively  expensive. 
The  basic  idea  of  domain  decomposition  is  to  solve  the  system  (2.3)  by  preconditioned  conjugate 
gradient  (PCG)  methods.  In  applying  the  PCG  method,  one  needs  to  evaluate  the  matrix-vector 
product  Cw  for  a  given  vector  w.  From  (2.4),  it  is  easily  seen  that  each  evaluation  of  Cw  requires 
the  solution  of  two  subdomain  problems.  For  example,  the  product  -  AfgAjfj1  A13W  can  be  computed 
by  solving  the  discretized  version  of  (2.1)  on  Hi  with  homogeneous  right  hand  side  (i.e.,  the  Laplace 
equation)  and  the  boundary  condition  u  =  w  on  O3,  and  then  taking  the  solution  on  the  first  row 
of  grid  points  above  Q3. 

3.  The  Eigen  Decomposition  of  the  Capacitance  Matrix 

In  order  to  understand  the  performance  of  various  preconditioners  for  C,  it  is  necessary  to 
first  analyze  the  eigen-structure  of  C.  It  turns  out  that  for  the  model  problem,  an  exact  eigen- 
decomposition  of  C  can  be  derived  by  the  use  of  Fourier  analysis. 

Define  the  vectors  Wj,j  =  1, 2,  •  •  • ,  n  by 


wj  =  \/2h(sin  jnh,  sin  2jich,  ■  ■  •  ,sin  njnh)T . 

(3.1) 

Consider  the  product 

Cwj  =  A33IUJ  A ^3 A j ^  A 1 3 Wj  A23A22  A23U7j*  . 

(3.2) 

Let  us  consider  the  term  -Af3AJ'11Ai3tyJ  first.  As  mentioned  earlier,  this  requires  the  solution  of 
the  discrete  Laplace  equation 

AhV  =  0  on  fit  (3.3) 

with  boundary  condition 

v  Wj  on  n3 

and  v  =  0  on  9nj/fl3  . 

Consider  a  solution  vector  v(x,  y)  of  the  form 

(3.4) 

v(ih,kh)  =  dk(wj)i  =  dkV2hsintjnh  , 

(3.5) 

where  0  <  i  <  n  +  1  and  0  <  k  <  mi  4-  1. 
The  boundary  condition  (3.4)  implies  that 


do  =  1  and  dm,+i  =  0. 


Substituting  (3.5)  into  (3.3),  we  get 

(dt-i  -  2 dk  +  d*+i)sin  ijith  +  d*(sin(i  -  1  )jirh  -  2sin  ijirh  +  sin(«  +  1  )jrh)  =  0. 


(3.6) 


It  follows  that  the  d*’s  satisfy  the  following  difference  equation  : 

dt-i  -  (2  +  <fj)dk  +  dk+i  —  0 


with  the  boundary  condition  (3.6)  and  where 


_  .  .  j J*h 
Os  s  4  sin  — — . 
'  2 


The  roots  of  the  characteristic  polynomial  corresponding  to  (3.7)  are 


r+  =  l  +  Zf  +  \/'i  +  -f 


and  r-  =  1  +  y  -  Joj  + 


The  general  solution  to  (3.7)  is  therefore  given  by 


d*  =  Cjr*+C2ri 

The  constants  c\  and  c?  can  be  found  by  imposing  the  boundary  conditions  (3.6)  which  give 


ci  ■  - 


r«*t+i 


_*»!+!  _  ,»|+I 
'+  r- 


and  cj  = 


r"*i+i 


_"»i+i  mi+i 
r+  ~  r- 


We  therefore  have 


- ^4j3i4n*v4j3Wy  s  d\Wj 


where 


By  a  similar  computation,  we  have 


r_  -  r+7?*,+I 

-<-nr&rH  , 


r_ 

S  ^ 


.r,-ia _ 


~^23^22  ^23  wj  “  (“ 


i -nr" 


Finally,  it  can  easily  be  verified  that 


{Az3Wj)i  =  sin(i  -  1)jjt/»  -  4  sin  +  sin(t  +  l)jwh, 


and  therefore 


AuWj  =  (-2  -  Oj)wj. 


Combining  (3.9),  (3.11)  and  (3.12),  we  obtain 


where 


Cwj  —  XjWj 


r_  —  r+'jfj*l+1  r_  —  r+7^**+1 
Ay  m  (  -  2  -  <Tj  +  (  .  _  m.4-1  )  +  (- 7~^TI  ))» 

1  < *  1  T» 


which  after  some  simplification  gives 


l  +  ^m1+l  1  +  I 

i  (j  _  ~«i+i  +  i  _  V  ^  +  4‘ 


(3.13) 


(3.14) 


Therefore  (Ay,  uiy)  is  an  eigenpair  of  C.  Since  the  vectors  tcy,  j  —  1, 2,  •  •  • ,  n  are  orthonormal, 
we  have  the  exact  eigen-decomposition  of  C. 

Theorem  3.1.  For  the  problem  (2.1),  the  capacitance  matrix  C  can  be  decomposed  as 

C  =  W\WT 

where 

A  «  diag(X1,X2y  -,X„), 

and 

W  =  {wi,W2,  • 

Furthermore,  W  is  orthogonal. 

Note  that  the  products  Cw  and  C~tw  can  be  computed  by  the  Fast  Fourier  TYansform  (FFT) 
(specifically,  the  sine  transform).  It  follows  that  for  the  model  problem,  there  is  no  need  to  use 
a  preconditioned  conjugate  gradient  method  for  the  interface  system:  we  can  solve  it  directly  by 
FFT. 

4.  Comparison  of  Preconditioners 

Among  the  preconditioners  proposed  so  far  for  C,  two  typical  ones  are  : 

MD  =  Wdia9(X?,X°,-  .  ,xZ)WT  (Dryja  [5])  (4.1) 

where 

A f  s  -2^37  (4.2) 

and 

Mg  =  Wdiag{ Af ,  A?,  •  •  • ,  A %)WT  (Golub-Mayers  [8])  (4.3) 

where 


Since  we  know  the  exact  eigen-decomposition  of  C  for  the  model  problem,  we  can  compare  the 
performance  of  these  two  preconditioners.  Specifically,  we  are  interested  in  the  spectral  condition 
numbers  K(MplC)  and  K(M^lC),  since  these  play  a  major  role  in  the  convergence  rate  of  the 
PCG  method. 


It  can  easily  be  verified  that 


r+  >  1,  0  <  r.  <  1  and  7,-  <  1. 


Therefore,  it  follows  immediately  from  (3.14)  that 


s  ~2\p  + V 


mj— 00 
ma—oo 


In  other  words,  for  fixed  h,  as  the  “aspect  ratios”  of  fti  and  Q2  tend  to  infinity,  A?  — ►  A j. 
that 

lim  K(Mq1C)  =  1. 


mi— 00 
m2— 00 


Next  let  us  consider  the  case  where  It  and  {2  are  fixed  but  h  —*  0  (i.e.,  n,  mi  and  m2 
the  eigenvalues  of  A I£lC  be  denoted  by  ftj,  where  by  (4.4)  and  (3.14) 

1,1  +  7“,+1  ,  l  +  7p+\ 

~  2  l  -  -vmi+1  l  -  V**+1 
*  'j  4  1i 

It  can  be  verified  that  (ij  is  a  decreasing  function  of  /,  and  therefore 

K(MclC)  =  !±  . 

/*n 

Consider  the  term  7™,+1  in  the  limit  as  h  — ►  0  for  It  and  I2  fixed.  We  have 


It  follows 


00).  Let 


where 


Making  use  of  the  fact  that 


we  have 


where 


Therefore,  we  have 


7r+i  -  (i  - 


6s  yfcH^j/4 

l  +  frj/2 


lim(l  +  xf(x))l/*  =  elim«— 0  /(*), 

x—O 


lim  7P,+l  =  e_2Ml 
h-o  1 


60  =  lim  7- - 777 

h-o  (1  +  6)h 

2sinTr 

=  lim  - ; — *~ 

h—0  h 


,.  l,l  +  e"2*'«  ,  l  +  e-2W», 

i-oMl  ~  2(l-e-2">  +  1-e-2"^  ' 


(4.10) 


(4.11) 


On  the  other  hand,  it  is  easy  to  verify  that 


lim  on  —  4 

A— 0 


and 

Therefore, 
and  so 


lim 
A— o 


3-2y/2 
3  +  2>/2 


<  1. 


lim'rr,+1  -  0, 

A— o 


and  lim  j”i+1  =  0 
A— o  " 


lim  un  =  1. 
0 


Combining  (4.11)  and  (4.12)  we  have 


Theorem  4.1.  For  any  lj  and  I2, 


(4.12) 


lim  *(MC->C)  "  +  TTT^) 


#-! 


1.1  +  e_2,i*  .  1  +  e~2*>2 , 


Note  that  the  limiting  value  of  K(Mq1C)  is  independent  of  fc.  Theorem  4.1  also  shows  that 
in  the  limit  h  — *  0,  the  rate  with  which  K(MpC)  — ►  1  is  exponential  in  l\  and  I2.  On  the  other 
hand,  if  1 1  and  h  are  small  then  the  limiting  value  of  K[Mq1C)  is  given  by 


S3  -£<£  +  £>.• 

‘fh—° 


and  thus  grows  like  £  and  jk. 

Next,  we  consider  K(MpC).  The  eigenvalues  Hj  of  MpC  are  given  by 


1  +  ^l+l 
1  _  ^7*1+1 


+ 


1  +  "Ip*1 
1  -  ip*1 


Moreover,  we  have 

The  following  result  follows  immediately. 
Theorem  4.2.  For  any  h,  1 1  and  h 


and 


1  K(Mp  C)  <  ^ 
>/2  "  K(MpC)  ~ 


lim  K(MpC)  =  v/2  . 


(4.13) 


We  thus  see  that  the  limiting  value  of  K(Mq1C)  is  also  independent  of  h,  but  this  value  does 
not  tend  to  unity  even  as  l\  and  h  tend  to  infinity.  This  is  a  fundamental  limitation  of  Mp. 

Plots  of  K{M£1C)  and  K(MplC)  for  various  values  of  h  and  l  in  the  case  where  li  =  =  l  are 

given  in  Figures  1  and  2.  It  is  seen  that  when  the  aspect  ratio  is  approximately  one  or  larger,  Mg 
is  a  pretty  good  preconditioner  for  C.  But  for  small  aspect  ratios,  both  Mg  and  Mp  deteriorate 
rapidly  as  preconditioners  for  C.  It  is  also  interesting  to  note  that  K  (Mq1C)  tends  to  its  asymptotic 
value  even  for  very  large  values  of  h  and  that  K(Mq1C)  has  a  smaller  value  for  larger  values  of  h. 

Bjorstad  and  Widlund[l]  proposed  a  preconditioner  Mb  (referred  to  as  the  “excellent  method” 
in  their  paper)  which  corresponds  to  the  capacitance  matrix  C  assuming  l\  —  h-  By  employing 
special  symmetry  in  the  case  where  lj  =  I2,  they  derived  a  method  for  computing  CMglw  exactly 
for  a  given  vector  w  by  solving  Neumann  problems  on  the  subdomains.  Obviously,  this  method 
is  exact  when  /j  —  /2  and  in  general  it  should  be  slightly  better  the  Golub-Mayers  preconditioner 
because  it  takes  into  account  the  aspect  ratio  of  one  of  the  subdomains. 

5.  Concluding  Remarks 

Our  analysis  provides  some  insight  into  the  structure  of  the  capacitance  matrix  system,  which 
is  central  to  the  domain  decomposition  method.  The  exact  eigen-decomposition  of  this  matrix  in 
the  simple  case  considered  in  this  paper  allows  us  to  compare  the  performance  of  various  precon¬ 
ditioned.  It  also  illustrates  somewhat  more  clearly  the  origin  of  the  pseudo  differential  operator 
used  in  Dryja’s  original  precondition?:  and  the  way  in  which  the  Golub-Mayers  preconditioner  is 
an  improvement.  Our  analysis  reveals  for  the  first  time  the  dependence  of  the  performance  of  these 
preconditioners  on  the  aspect  ratios  of  the  subdomains.  Based  on  the  results  here,  it  is  straightfor¬ 
ward  to  derive  a  preconditioner  for  irregular  domains  which  takes  into  account  che  aspect  ratios  of 
the  subdomains.  Furthermore,  the  knowledge  of  the  eigen-decomposition  of  C  makes  it  possible  to 
construct  a  direct  domain-decomposed  fast  Poisson  solver  on  a  rectangle  divided  into  many  strips 
and  boxes.  We  referred  the  interested  readers  to  [4j. 


Acknowledgement:  The  author  thanks  Diana  Resasco  for  helpful  discussions  and  producing 
the  plots  in  this  paper. 


Figure  1:  Dependence  of  K{Mq1C)  on  h  and  aspect  ratio  l 
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Figure  2:  Dependence  of  K(Af£lC)  on  h  and  aspect  ratio  l 
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